{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Slightly compressible flow\n",
    "In this tutorial we investigate how to use PorePy to solve slighlty compresible flow. <br>\n",
    "\n",
    "\n",
    "**Note**: This tutorial will focus on a 2d domain, however most of the code works for 1d, 2d, and 3d domains.<br>\n",
    "\n",
    "Let is $\\Omega$ a regular domain with boundary $\\partial \\Omega$. The boundary can be divided in two non-overlapping parts useful to impose Dirichlet ($\\partial \\Omega_d$) and Neumann ($\\partial \\Omega_n$) boundary conditions. We indicate with $\\mathbf{n}$ the outward unit normal vector of $\\partial \\Omega$.<br>\n",
    "The single-phase flow can be written in primal formulation as:\n",
    "\n",
    "$$ \\phi\\frac{\\partial\\rho}{\\partial t}- \\nabla \\cdot \\rho K \\nabla p = \\rho f $$\n",
    "with boundary conditions on $\\partial \\Omega_d$ and $\\partial \\Omega_n$:\n",
    "$$ p = p_b \\qquad - K \\nabla p \\cdot \\mathbf{n} = u_b$$\n",
    "\n",
    "Where $\\phi$ is the porosity, $\\rho$ is the fluid density, $f$ is a scalar source/sink term, $K$ is the permeability matrix, $p_b$ is the pressure at the boundary (Dirichlet condition), and $u_b$ is the flux at the boundary (Neumann condition).<br>\n",
    "\n",
    "As a relationship between pressure and density we use $c_p\\rho = \\text{d}\\rho/\\text{d}p$. Assuming slightly compressible flow (e.g., $\\nabla\\rho\\cdot K\\nabla p \\ll 1$) we can write conservation of mass as\n",
    "\n",
    "$$c_p\\phi\\frac{\\partial p}{\\partial t} - \\nabla \\cdot K \\nabla p = f $$\n",
    "\n",
    "We now discretize in time using backward Euler and time step $k$:\n",
    "$$c_p\\phi\\frac{p^{k+1} - p^k}{k} - \\nabla \\cdot K \\nabla p^{k+1} = f^{k+1} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import modules\n",
    "\n",
    "First we have to import the modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import porepy as pp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define problem setup:\n",
    "We define our problem in the unit square. We let our data inherit from the compressible data assigner base class to set default values. We set different parameters in the matrix and fractures. For the matrix we use use default parameters for everything except the initial pressure. Note that we let the Fracture data inherit from the matrix data. In this way, the initial pressure (or any other parameters you may set for the matrix) is equal the parameter for the matrix unless it is overloaded in the fracture data class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MatrixDomain(pp.SlightlyCompressibleDataAssigner):\n",
    "    def initial_pressure(self):\n",
    "        return np.ones(self.grid().num_cells)\n",
    "\n",
    "class FractureDomain(MatrixDomain):\n",
    "    def permeability(self):\n",
    "        kxx = 1000 * np.ones(self.grid().num_cells)\n",
    "        return pp.SecondOrderTensor(2, kxx)\n",
    "    \n",
    "    def aperture(self):\n",
    "        return np.power(0.001, 2 - self.grid().dim)\n",
    "\n",
    "class IntersectionDomain(FractureDomain):\n",
    "    def source(self, t):\n",
    "        assert self.grid().num_cells == 1, 'Assumes Intersection domain only has 1 cell'\n",
    "        f = .4 * self.grid().cell_volumes  # m**3/s\n",
    "        return f * (t < .05)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to create the grid bucket and assign data to it. The important thing to notice is that the data classes should be stored with the keyword 'flow_data' (or more general, solver.physics + '_data')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate grid bucket\n",
    "nx = 12\n",
    "ny = 12\n",
    "frac1 = np.array([[0, 1], [.5, .5]])\n",
    "frac2 = np.array([[.5, .5], [0, 1]])\n",
    "fracs = [frac1, frac2]\n",
    "physdims = [1, 1]\n",
    "gb = pp.meshing.cart_grid(fracs, [nx, ny], physdims=physdims)\n",
    "\n",
    "# Initate data\n",
    "for g, d in gb:\n",
    "    if g.dim == 2:\n",
    "        d['flow_data'] = MatrixDomain(g, d)\n",
    "    elif g.dim == 1:\n",
    "        d['flow_data'] = FractureDomain(g, d)\n",
    "    elif g.dim == 0:\n",
    "        d['flow_data'] = IntersectionDomain(g, d)\n",
    "    else:\n",
    "        raise ValueError('Unkown grid-dimension %d' %g.dim)\n",
    "\n",
    "# We loop over the edges and assign the coupling permeability\n",
    "for e, d in gb.edges():\n",
    "    d['kn'] = 1000\n",
    "\n",
    "# Set time stepping parameters:\n",
    "dt = .005\n",
    "T = 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now initialize the problem and solve it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Time stepping using 20 steps\n",
      "Step 0 out of 20\n",
      "Step 1 out of 20\n",
      "Step 2 out of 20\n",
      "Step 3 out of 20\n",
      "Step 4 out of 20\n",
      "Step 5 out of 20\n",
      "Step 6 out of 20\n",
      "Step 7 out of 20\n",
      "Step 8 out of 20\n",
      "Step 9 out of 20\n",
      "Step 10 out of 20\n",
      "Step 11 out of 20\n",
      "Step 12 out of 20\n",
      "Step 13 out of 20\n",
      "Step 14 out of 20\n",
      "Step 15 out of 20\n",
      "Step 16 out of 20\n",
      "Step 17 out of 20\n",
      "Step 18 out of 20\n",
      "Step 19 out of 20\n"
     ]
    }
   ],
   "source": [
    "solver = pp.SlightlyCompressibleModel(gb, end_time=T, time_step=dt)\n",
    "solution = solver.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can plot the pressure field at the end time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVwAAADxCAYAAACH4w+oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8lNW9+PHPdyYLBDDsLoQ9aBta\nijRsYntZWkDUoK8ignWrSLVCiz+tqPVCU1QEraiXpRalIlgDVF+aVFlcuLGvctki3IrgFaKyJEgh\nIFtolpn5/v6YSZoMIXmGTCaT4ft+vZ4XM8+cOec884Rvnpx5zvmKqmKMMabhuRq7A8YYc6GwgGuM\nMRFiAdcYYyLEAq4xxkSIBVxjjIkQC7jGGBMhFnCNMRcEERktIp+LSL6IPFLD64kisjLw+mYR6Rb0\nehcROS0iv3ZaZzALuMaYmCcibmAhcA2QBkwUkbSgYpOAb1Q1FXgOmBv0+jxgTYh1VmMB1xhzIRgA\n5Kvql6paBqwAxgaVGQu8Gnj8BjBCRARARG4AvgJ2hlhnNXEhdtqmpRljnJL6vDlVRM84LPu1PxCW\nVNm1WFUXV3neCThQ5XkBMDComsoyquoRkRNAOxEpAR4Gfgz8uqbytdRZTagB1xhjIuIMcI/DsplQ\noqrpDdSVTOA5VT0duOA9bxZwjTFRSQhrgCoEOld5nhLYV1OZAhGJA5KBo/ivWseJyNNAa8AXuOr9\n2EGd1VjANcZEJRfQPHzVbQV6iUh3/EFxAnBLUJkc4A5gIzAOWK/+1b1+UFFARDKB06q6IBCU66qz\nGgu4xpioJEB8mOoKjMlOBdYBbuBPqrpTRGYBeaqaAywBlotIPnAMfwANuc7a3iMhLs9oX5oZY5yq\n14BnVxGt88bWgPvg4wYcww0bu8I1xkSlcF7hRgsLuMaYqBTmL82iQqwdjzEmRtgVrjHGREiY71KI\nChZwjTFRya5wjTEmgmItQMXa8RhjYoRd4RpjTITYXQrGGBMh9qWZMcZEiA0pGGNMhNiQgjHGRIhd\n4RpjTITYFa4xxkSIXeEaY0yECHaXgolBHo+HuLjG/VGIhj6Y6CJAvNMfCU9D9iR8LE16DOvWrRtP\nPfUUaWlptGnThp/97GeUlJSQm5tLSkoKc+fO5ZJLLuFnP/sZAO+88w59+/aldevWXHXVVXzyySeV\ndc2dO5dOnTrRqlUrrrjiCj788EMAtmzZQnp6OhdddBEXX3wxDzzwAEBlG8H9+eCDDwDIzMxk3Lhx\n3HrrrVx00UUsXboUn8/HnDlz6NmzJ+3atWP8+PEcO3YsEh+ViUIiEBfnbGsqLODGuD//+c+sW7eO\nL774gt27d/PEE08AcOjQIY4dO8a+fftYvHgx27dv56677uKPf/wjR48e5Z577iEjI4PS0lI+//xz\nFixYwNatWzl16hTr1q2jW7duAEybNo1p06Zx8uRJvvjiC8aPH++4b9nZ2YwbN47jx4/z05/+lPnz\n5/P222/z0UcfcfDgQdq0acOUKVMa4mMxTYAIxLudbU2FBdwYN3XqVDp37kzbtm157LHHyMrKAsDl\ncvG73/2OxMREmjdvzuLFi7nnnnsYOHAgbrebO+64g8TERDZt2oTb7aa0tJRdu3ZRXl5Ot27d6Nmz\nJwDx8fHk5+fz/vvvc9VVV9G3b1+Ki4u58847KS8vr7VvgwcP5oYbbsDlctG8eXNefPFFnnzySVJS\nUkhMTCQzM5M33ngDj6eJ/L1owircV7giMlpEPheRfBE5K3uPiCSKyMrA65tFpFtg/wAR+d/A9g8R\nubHKe/aKyI7Aa3l19cECbozr3PnfWZy7du3KwYMHAejQoQPNmjWrfG3fvn08++yztG7dunI7cOAA\nBw8eJDU1leeff57MzEw6duzIhAkTKutZsmQJu3fvZuLEiRQVFTFx4kSmT5/Oj3/8Y+Lja/+OuWrf\nKvpw4403Vrb/7W9/G7fbzT//+c9wfRymCRGB+ERnW911iRtYCFwDpAETRSQtqNgk4BtVTQWeA+YG\n9n8KpKtqX2A08MdAxt4Kw1S1r5OcahZwY9yBAwcqH+/fv5/LLrsMAJHq+f06d+7MY489xvHjxyu3\nM2fOMHHiRABuueUW/v73v7Nv3z5EhIcffhiAXr16kZWVxeHDh5k3bx7Z2dls3ryZSZMmcebMmcr6\nvV4vR44cqdZmTX1Ys2ZNtT6UlJTQqVOn8H0gpumouBHXyVa3AUC+qn6pqmXACmBsUJmxwKuBx28A\nI0REVPWMqlb8mdWMeiTTtYAb4xYuXEhBQQHHjh3jySef5Oabb66x3OTJk3nxxRfZvHkzqkpxcTHv\nvvsup06d4vPPP2f9+vWUlpbSrFkzmjdvjsvl/9F57bXXOHLkCC6XqzKAnj59mi5dulBSUsK7775L\neXk5TzzxBKWlpbX29d577+Wxxx5j3759ABw5coTs7OwwfhqmSQlvwO0EHKjyvCCwr8YygQB7AmgH\nICIDRWQnsAO4t0oAVuA9EflYRH5eVycs4Ma4W265hZEjR9KjRw969uzJf/7nf9ZYLj09nZdeeomp\nU6fSpk0bUlNTWbp0KQClpaU88sgjtG/fnksuuYTDhw/z1FNPAbB27Vp69+5Ny5Ytufvuu7n//vu5\n7bbbmD17NosWLeLuu++mU6dOtGjR4qy7FoJNmzaNjIwMRo4cSatWrRg0aBCbN28O6+dhmhjnAbe9\niORV2eoMfqFQ1c2q2hvoDzwqIhXjcVeraj/8QxVTROSHtdUjqiFdHZ/3pbSJvG7duvHyyy/zox/9\nqMHbWrZsGdnZ2bz55pt4vV6uuuoqnnrqKYYPH97gbZuoJXUXObf0ZqJ5XRw2tIePaxtDFZHBQKaq\njgo8fxRAVZ+qUmZdoMzGwBjtIaCDBgVJEVkPTFfVvKD9mcBpVf39ufrRhO5gM9Hs9ttv5/bbbwfA\n7Xbblampv/AuprAV6CUi3YFCYAJwS1CZHOAOYCMwDlivqhp4zwFV9YhIV+BbwF4RaQG4VPVU4PFI\nYFZtnbCAa4yJTgI4uAPBiUCwnAqsA9zAn1R1p4jMAvJUNQdYAiwXkXzgGP6gDHA18IiIlAM+4D5V\nLRKRHsBbge8u4oDXVXVtrYdkQwrGmAZSvyGFFqJ5wTdunauhvNqHFKKFXeEaY6JTDK7PGGOHY4yJ\nKU1o2q4TFnCNMdHJrnCNMSZCLOAaY0yEhPEuhWhhAdcYE51i8ArXpvY6tHbtWq644gpSU1OZM2fO\nWa+XlpZy8803k5qaysCBA9m7d2+DtTVv3jzS0tLo06cPI0aMqFx7oCHaqvDmm28iIuTl1bkCXb3b\nW7VqFWlpafTu3Ztbbgm+Nz18be3fv59hw4Zx5ZVX0qdPH1avXn3ebd1111107NiR73znOzW+rqr8\n6le/IjU1lT59+rBt27bzbuuCEd61FKKDqoayXZA8Ho/26NFDv/jiCy0tLdU+ffrozp07q5VZuHCh\n3nPPPaqqmpWVpePHj2+wttavX6/FxcWqqrpo0aIGbUtV9eTJk/qDH/xABw4cqFu3bj2vtpy2t3v3\nbu3bt68eO3ZMVVX/+c9/NlhbkydP1kWLFqmq6s6dO7Vr167n1Zaq6kcffaQff/yx9u7du8bX3333\nXR09erT6fD7duHGjDhgw4LzbakJCjS/Vtu+3QfUmZxv+yQv1ai8Sm13hBikvL8fj8aBVJoRs2bKF\n1NRUevToQUJCAhMmTDhrFavs7GzuuOMOAMaNG8eHH35YrQ6nnLQ1bNgwkpKSABg0aBAFBQUht+O0\nLYAZM2bw8MMPV1s/t6Hae+mll5gyZQpt2rQBoGPHjg3Wlohw8uRJAE6cOFG5dOX5+OEPf0jbtm3P\n+Xp2dja33347IsKgQYM4fvw4X3/99Xm3d0GIwStcC7hBvF4vX331VbWgW1hYWG2x7JSUFAoLC6u9\nr2qZuLg4kpOTOXr0aMjtO2mrqiVLlnDNNdeE3I7TtrZt28aBAwe49tprz6uNUNvbvXs3u3fvZsiQ\nIQwaNIi1a2udKVmvtjIzM3nttddISUlhzJgxzJ8//7zaCld/TJCKL82cbE2EBdwgPp+PgoICPB4P\nX3zxxXldpUbKa6+9Rl5eHg899FCD1O/z+XjggQd49tlnG6T+mng8Hvbs2UNubi5ZWVlMnjyZ48eP\nN0hbWVlZ3HnnnRQUFLB69Wpuu+02fD5fg7RlzoNd4ca+o0ePUlxcTHl5OQcOHCA/P59LL720WuaE\ngoKCs7IQdOrUqbKMx+PhxIkTtGvXLuT2q9ZzrrYAPvjgA5588klycnJITDy/X/F1tXXq1Ck+/fRT\nhg4dSrdu3di0aRMZGRnn/cWZk2NLSUkhIyOD+Ph4unfvzuWXX86ePXsapK0lS5ZUJr0cPHgwJSUl\nFBUVhdxWuPpjgljAjX0dOnQgMTGRvLw8vF4vBQUF9OnTh/bt2/PVV19RVlbGihUryMjIqPa+jIwM\nXn3Vn53jjTfeYPjw4WelkHGif//+7Nmzp9a2tm/fzj333ENOTs55j3E6aSs5OZmioiL27t3L3r17\nGTRoEDk5OaSnn98aIU6O7YYbbiA3NxeAoqIidu/eTY8ePRqkrS5dulSme//ss88oKSmhQ4cO53Vs\ndcnIyGDZsmWoKps2bSI5OZlLL720QdqKKTEWcJtQVyMnLi6Ovn37smHDBuLj44mLi+Omm25i5MiR\n+Hw+7rrrLnr37k3Pnj157rnnyMjIYNKkSdx2222kpqbStm1bVqxYcd5tL1iwgFGjRuH1eivbmjlz\nJu+88w7btm3joYce4vTp09x0002AP3Dk5OQ0SFvh5KS9UaNG8d5775GWlobb7eaZZ545r78Uamsr\nPT2dRYsWMW/ePCZPnsxzzz2HiLB06dLz+iUJMHHiRHJzcykqKiIlJYXf/e53lVmL3377bdasWcPq\n1atJTU0lKSmJV1555bzauaAIMbeWgi3PGKSkpISNGzdy1VVXsWHDBkpKSujYsSMnT56kc+fOdOnS\nhYSEBESE9PT0et+X6pTP52PAgAERay+Wj83r9TJw4MCItBfpY4sy9Vue8RLRvFsdNvSsLc/Y5IkI\nzZs3Jzk5mUOHDrF//34ApkyZQlFREUePHqVfv34R6cugQYPYu3dvxNqL5LH169ePffv2Ray9yZMn\ns2vXroi0993vfpf9+/dXtrV///4GGyeOOTa198LUpUsX9u3bx5kzZ/jqq6/4y1/+gsfj4eabb2bl\nypUR6UNhYSH79u1j2bJlEWkvksd26NAhDh48yPLlyyPS3vbt2xk1ahQvv/xyg7d1+PBhjhw5wvLl\nyxGRiOSXixkxOLU3xg6n4bjdblq0aEFJSQmfffYZPp+PjVs20759e8d1SLwbLfc6KuuKd+Oroey5\n2jtX+ZrLuvCV1337U0Vb7ngXXgflAUbPGsDamVsclQ3uS12fZSjH6KT822+/Xfk4lHMTavmki1ry\n+eef8+CDD/J///d/pKen0759+/O+x/iCYQH3wlYxxNCvXz88Hg9nTp1mpJ49M+tc3pOxXK+rHJX9\nq4znNl3suO7l8nN+qU87KjtfpjNDf+O47sdlNk/rLx2VdeNl+IyBjuueLvMd9+Vxme34GMF/nE4/\nw+Xyc8fnBvznx+m5f0/G0rdvXz788EOuvvrqC3U8N3QxGHDttjBjGpktfFMLt8PNAREZLSKfi0i+\niDxSw+uJIrIy8PpmEekW2D9ARP43sP1DRG50WmcwC7jGNLI777yz1uGFNWvWsGfPHvbs2cPixYv5\nxS9+EcHeNaIwTnwQETewELgGSAMmikhwispJwDeqmgo8B8wN7P8USFfVvsBo4I8iEuewzmos4BrT\nyGzhm3MI71oKA4B8Vf1SVcuAFcDYoDJjgVcDj98ARoiIqOoZVfUE9jfj37fHOqmzGgu4xkS5C3bh\nm9CucNuLSF6V7edBtXUCDlR5XhDYV2OZQIA9AbQDEJGBIrIT2AHcG3jdSZ3VxNiQtDEmZoT2pVlR\nQ058UNXNQG8R+TbwqoisOZ967ArXmCh3wS58UzG1NzxfmhUCnas8Twnsq7GMiMQByUC1NVZV9TPg\nNPAdh3VWYwHXmCh3wS58E97VwrYCvUSku4gkABOA4AVIcoA7Ao/HAetVVQPviQMQka7At4C9Duus\nxoYU6kHi3bwntY6RVy8f5+avMt5hWRfLzxqGqr38fJnuqKwrzsXjMttx3a44Ybo4W5z7vtyxLBrq\n/N7kUPoSyjFWlHf6GYZybirKOz33El/7JVhtC9/ce++9jBkz5sJc+Ebwf0UVBqrqEZGpwDr818R/\nUtWdIjILf3qeHGAJsFxE8oFj+AMowNXAIyJSDviA+1S1CKCmOmvrhwXcetByb8g3y4/XpY7KrpI7\nmaQLHNe9RKbyoD7uqOyzMoNMfdhx3ZkyN6SJD07Lgn/ig9O+ZMpcx8cI/uN0+hkukamOzw34z08o\nk1hqk5WVVevrIsLChQsd9y1mhHm1MFVdDawO2jezyuMS4KYa3rccqHHeeU111sYCrjEmOsXgTLMY\nOxxjTEyJsQgVY4djjIkZMbgAuQVcY0x0siEFY4yJEFuA3BhjIsSucI0xJkIs4Ma+0tJSfD5n2Q2M\nMQ0oBgOuZe0NcuTIETZv3kyHDh04deoUbrebq666iv/5n/+pLDNkyBA8Hg8tW18UWlqWOBfqcRbM\nQykL/hlbPoflQynrLy/4PM5O/ZTc61k49K8h1N1w/W7Iz1vi3KjH2bmXeDenj58E4Oqrr76QFhCv\nX9be3qJ5DtPqyXcta2+T1KpVK1q0aEHXrl35+OOPAX8G25pouTfkNDihzHwKNZ3Mb3SGo7Kz5XGe\n0vsd1/2oPM+zep+jsm48jssCPCiLHPflUXne8TGC/zhDSTsU6sy+UNL3mPMQg1e4MXY44dOmTRuS\nkpLwer0UFhZSXFxMQkICcXH2kRkTEXaXwoXH7XbTp08f/v73v1NWVkZpaSkbN24kxKEYc4HbsmUL\n8fHxjd2NpsWucC9cLpeLZs2aoaoMGDAAr9f52K0x6enp3HjjjZYmPRQxGHBtPdwQiQhutxuXyz46\n45zL5SI7O5tvfetb5OXlnRVs165dyxVXXEFqaipz5sw56/379+9n2LBhXHnllfTp04fVqx0vUNV0\nhXc93KhgUcOYRub1epkyZQpr1qxh165dZGVlsWvXrmplnnjiCcaPH8/27dtZsWIF993n/IvJpkzd\nzramogn9bjAmNm3ZsoXU1FR69OgBwIQJE8jOziYt7d8Zt0WEkyf9t5adOHGCyy67rFH6GknqgrIw\nLUAeLSzgGtPIasrKu3nz5mplMjMzGTlyJPPnz6e4uJgPPvgg0t2MOBXwuJ3+Ed40JitZwK0HV7w7\n5DQ4S2Sq47KhpJNxxbmYLc6yIbjiXDwqz4dQt/CgLHJUdmrudSwY+k5IdTvtSyjHCKF9hqGcm4ry\nTs+9q44UO05kZWVx55138uCDD7Jx40Zuu+02Pv3005j+LkFF8Dq+DbOsQfsSLhZw68FX7g15ckIo\naXBCvck/lAkEoU5OWKCTHJX14HZcFmCqLHHcl1AmSUBoEyVmy+Mhp+8JZVJFbZxk5V2yZEnlF22D\nBw+mpKSEoqIiOnbs6LjPTZHXHb4BWhEZDbyAf5Xdl1V1TtDricAy4Pv4s/XerKp7ReTHwBwgAX9k\nf0hV1wfekwtcCvwrUM1IVT18rj7E7q9HY5qI/v37s2fPHr766ivKyspYsWIFGRkZ1cp06dKFDz/8\nEIDPPvuMkpISOnTo0BjdjRhF8OJ2tNVFRNzAQuAaIA2YKCJpQcUmAd+oairwHDA3sL8IuF5Vv4s/\nq29wfrOfqmrfwHbOYAsWcI1pdHFxcSxYsIBRo0bx7W9/m/Hjx9O7d29mzpxJTo4/6/azzz7LSy+9\nxPe+9z0mTpzI0qVLEanXUgVRTxE8uB1tDgwA8lX1S1UtA1YAwWmXxwKvBh6/AYwQEVHV7ap6MLB/\nJ9A8cDUcMhtSMCYKjBkzhjFjxlTbN2vWrMrHaWlpbNiwIdLdalSKUOZ8bm97Ecmr8nyxarXFLjoB\nB6o8LwAGBtVRWSaQVv0E0A7/FW6FnwDbVLW0yr5XRMQLvAk8obVMQ7WAa4yJShVDCg4VNfRqYSLS\nG/8ww8gqu3+qqoUi0gp/wL0N/zhwjWxIwRgTtcI1hgsUAp2rPE8J7KuxjIjEAcn4vzxDRFKAt4Db\nVfWLijeoamHg31PA6/iHLs7JAq4xJiqFeQx3K9BLRLqLSAIwAcgJKpOD/0sxgHHAelVVEWkNvAs8\noqqV4zoiEici7QOP44HrgE9r64QNKRhjopJ/SCE8ISowJjsVWIf/trA/qepOEZkF5KlqDrAEWC4i\n+cAx/EEZYCqQCswUkZmBfSOBYmBdINi6gQ+Al2rrh2V8CFJSUsLGjRurZXk4V8aHVq0vwhdCxoeG\nzcoQSt3OMziEWv5XudfyX0PfbZC6Q+93dGSTcMW7OWUZH0KWlt5c/5zX3VHZfvKZZXyIdb5yLzP0\nN47LPy6zydSHHZXNlLkNlpUhlIkM4J+cEEr5hqo7lEkSEHo2CafnBvznx+m5f1xmO67X/JuC0+GC\nJsMCrjEmSoVvSCFaxNbRGGNiRoi3hTUJFnCNMVHLAq4xxkSAXeFeAE6ePElJSQkFBQV4vd6YXv7O\nmGimCKUxlrbXAm6Q5s2bExcXh9frpby8HK/Xy4YNGzhz5gwulwuXy8XevXvx+ZrGgscmOuzfv99+\neYcoFq9w7ScgSHx8PHFxcXTt2pVmzZrRokULhgwZQlJSEgkJCbhcLuLi4ohzvDCyMf4Vwe6+++7K\nrL2jR49u7C5FvXAuzxgtLOA6VJGtNz4+npSUlAsip5QJn8suu4y1a9eed9ZegFWrVpGWlkbv3r25\n5ZZbItHtRhfGqb1RwWaaBQltplkrfOUNNRss1JlmDTljy3n5abljeGGo8xTe0dLvBp3ZF+/i1PFT\nQM0zzbxeL5dffjnvv/8+KSkp9O/fn6ysrGpJJPfs2cP48eNZv349bdq04fDhw00h20O9Zpr1TG+t\nc/J+6KjsePmrzTSLdb5yH0/rLx2Xny7zHZefLvMbLA1OqDPHpsoS/qB3Oirrxe24LMAvZGmD9juU\n2XcNeS5r4yRr70svvcSUKVNo06YNQFMItvVmY7jGmLCrKWtvYWH1lQN3797N7t27GTJkCIMGDTpr\nSCIW+e9SSHC0NRV2hWtME+DxeNizZw+5ubkUFBTwwx/+kB07dtC6devG7lqDCedqYdHCrnCNaWRO\nsvampKSQkZFBfHw83bt35/LLL2fPnj2R7mrE2V0KxpiwcpK194YbbiA3NxeAoqIidu/eXTnmG6ti\n8baw2LpeN6YJqpq11+v1ctddd1Vm7U1PTycjI4NRo0bx3nvvkZaWhtvt5plnnqFdu3aN3fUGFYtf\nmlnANSYK1JW1V0SYN28e8+bNq7WemTNn0rZtW+6/378O8GOPPUbHjh2ZNm1a+DvdwGJxaq8NKZiw\nOV1UwgfPfMLpopLG7soF66677mLZMn/SWJ/Px4oVK7j11lsbuVfnx4YUTDWjZw3AjfMUO/fljnVc\nfkpuBm48juuemnud4xk3v8y9znG94J/M4OSHesMrX/D29K0A/OihPo7rdupXudeGNKtoau51jj/D\nKbnXN9i5HD2r1kSuYdWtWzfatWvH9u3b+ec//8mVV17ZpIcewhlMRWQ08AL+/GMvq+qcoNcT8ac4\n/z7+bL03q+peEfkxMAdIAMqAh1R1feA93weWAs2B1cA0rWU2mV3h1sPamVsc/wb24mbR0GzHZRcO\nzcFLnONtwdB3HJeeP/SdkI7zhaGrHfV6yM96csPT/Rnys56OP5VQZqX919B3Q/hEvCwY+o7j0guH\n/rXBzuXamVtC/dGql7vvvpulS5fyyiuvcNddd0W07XAKZ9ZeEXEDC4FrgDRgooikBRWbBHyjqqnA\nc8DcwP4i4HpV/S7+rL7Lq7znD8BkoFdgq3WRDAu4Jmxatm/Gjx7qQ8v2zRq7Kxe0G2+8kbVr17J1\n61ZGjRrV2N05bxX34TrZHBgA5Kvql6paBqwAxgaVGQu8Gnj8BjBCRERVt6vqwcD+nUBzEUkUkUuB\ni1R1U+CqdhlwQ22dsCEFY2JMQkICw4YNo3Xr1rjdTWd8syYhDCm0F5G8Ks8Xq+riKs87AQeqPC8A\nBgbVUVkmkFb9BNAO/xVuhZ8A21S1VEQ6BeqpWmf1G6iDWMA1Jsb4fD42bdrEX/7yl8buSr0oQpnz\nabtFDb14jYj0xj/MMPJ867AhBWNiyK5du0hNTWXEiBH06tWrsbtTL+EcwwUKgc5VnqcE9tVYRkTi\ngGT8X54hIinAW8DtqvpFlfIpddRZjV3hGhND0tLS+PLLLxu7G2ER5rUUtgK9RKQ7/qA4AQheVDgH\n/5diG4FxwHpVVRFpDbwLPKKqGyr7p/q1iJwUkUHAZuB2oNal4ewK1xgTtcJ1H66qeoCpwDrgM2CV\nqu4UkVkiUjGPegnQTkTygQeARwL7pwKpwEwR+d/AVrE+5n3Ay0A+8AWwprZ+2BWuMSYqhXtqr6qu\nxn+vbNV9M6s8LgFuquF9TwBPnKPOPOA7TvtgAdcYE5UqxnBjiaXYCVJXih1VpX///ni9Xtp1bBdi\nip1oSSdjKXYiWbc73kXR4aMADB8+/KwUOzGsXil22qd31evzHnVUdqn8wlLsNFWqSnFxMR6PB6/X\ny44dOyguLkZVERE++eQTVBVfuY8Z+hvH9T4us8nUhx2VzZS5PKX3O677UXk+pHQyoaaqCaV8Q9Ud\nSsoc8B+n08/wUXne8bkB//lxeu4fl9l8+umn/PrXv67M2tu+ffsLImtDfcTiamH2pVmQo0ePUlxc\nzJ49e/B6vbhcLrp06UJSUhItW7akRYsW9O/fn379+jV2V00T0q9fP9avX1+vrL0Ab775JiJCXl7e\nOcvEklhbvMYCbpB27drRsmVL+vbtS2JiIvHx8SQnJyNSr7+OjDknr9fLlClTWLNmDbt27SIrK4td\nu3adVe7UqVO88MILDBwYPEEqNoX5PtyoYAHXmEZWNWtvQkJCZdbeYDNmzODhhx+mWbMLY62KMK+l\nEBUs4BrTyJxk7d22bRsHDhyvvXYHAAAVeklEQVTg2muvjXT3Gk3F1F4nW1PRdH41GHOB8vl8PPDA\nAyxdurSxuxJRsXhbmAVcYxpZXVl7T506xaeffsrQoUMBOHToEBkZGeTk5JCeHvV3QtVLUxoucCK2\njsaYJqhq1t5OnTqxYsUKXn/99crXk5OTKSr69wqBQ4cO5fe//33MB9tYvC3MAm49uOJdPC6znZeP\nc5Epc+suGCj7qDwfQt3Cg7LIcdmpsiSkup2W/2XudSFllAil7lCOsaK8088wlHNTUd7puXfF1/5V\niZOsvRciC7imGl+5j1/q047Lz5fpPKiPOyr7rMzgNzrDcd2z5fGQbvIPdQKB08kJHtwhT3wIZcJG\nqJNBnH6Gs+Vxx+cG/OfH6bmfL9PrLFNX1t6qcnNzHbUbC2wM1xhjIsCHi7IYS5NuAdcYE7VsSMEY\nYyLAxnCNMSZCFBvDNcaYCAlrip2oEFtHY4yJGbE4pGBrKRhjopIilJLgaHNCREaLyOciki8ij9Tw\neqKIrAy8vllEugX2txOR/xaR0yKyIOg9uYE6g3Od1dwHy/hQXV0ZHwCGDBmCx+OhVeuL8JV7Hdft\ninPh8zjLEBFK2dDrbrjsBlNzr2NBiBMfGi4rQ8N93hLnQp3WHe/m1PGTAFx99dWW8cGhhPTv6sV5\nbzkqWyC9as34ICJuYDfwY6AAfxbfiaq6q0qZ+4A+qnqviEwAblTVm0WkBXAl/txl31HVqVXekwv8\nOpDbrE42pFAPvnIvt+lix+WXy8+ZpAvqLggskakhT6oI5Sb/hsom4cbToFkZQp0MEsrkBKfnBvzn\nx+m5Xy4/d1yvqS6MQwoDgHxV/RJARFYAY4GqCw+PBTIDj98AFoiIqGox8HcRSa1vJ2xIwRgTlSrG\ncB1mfGgvInlVtuDfcp2AA1WeFwT21VgmkFb9BNDOQVdfCQwnzJA6MhXYFa4xJiopgtfn+Aq3qJGS\nSP5UVQtFpBXwJnAbsOxchS3gGmOikvqE0pKwTe0tBDpXeZ4S2FdTmQIRiQOSgaO19lG1MPDvKRF5\nHf/QhQXccPB4PJSVlbF582ZC/LLRXODy8vKIj49v7G40KaqC1xO2MdytQC8R6Y4/sE4AbgkqkwPc\nAWwExgHrtZb/6IGg3FpVi0QkHrgO+KC2TtgYrgMHDx6kuLiY8vJyEhMT6devH3379m3sbpkmpG/f\nvsyYMaMyTfro0aOrvV5X1t558+aRlpZGnz59GDFiBPv27YtU1xuPgtfjdrTVWZV/THYqsA74DFil\nqjtFZJaIVKx/uQRoJyL5wANA5a1jIrIXmAfcKSIFIpIGJALrROQT4H/xB/KXauuHXeGeg9frpays\njLKyMk6dOkVSUlJl5t74+HjL4mtCEhcXR05ODldfffVZKc4rsva+//77pKSk0L9/fzIyMkhLS6ss\nc+WVV5KXl0dSUhJ/+MMfmD59OitXroz0YUSUquApD9/EB1VdDawO2jezyuMS4KZzvLfbOar9fih9\nsCvcIOXl5ZSWlrJx40ZUlRYtWnDFFVdYgDUNxknW3mHDhpGUlATAoEGDKCgoaIyuRpjg88Y52poK\nC7hBiouLEREGDx5MYmKiBVrT4Jxk7a1qyZIlXHPNNZHoWuNSwON2tjURTedXQ4S0bt2ahIQE3O66\nT6Ir3h3STe0S52KJTK27YKCsk0wBlX2JczFbnGUsaMj0PVNyM1g4NCekukNJg+P0GCG0zzCUc1NR\n3um5d8WHLyC89tpr5OXl8dFHH4WtzqjlEyiJrRAVW0cTYb5yL9frKsfl/yrjGa9LHZVdJXeGPPMp\nlPQ9mfqw47ozZS5P6y8dlXXjdVwWYLrMd9yXTJkbchqcUGb2OT034D8/Ts/9X2V8ra/XlbW3wgcf\nfMCTTz7JRx99RGJibGVCOCdPY3cgvGxIwZhGVjVrb1lZGStWrDgrceT27du55557yMnJoWPHWtdH\niR3+BXGdbU2EXeEa08icZO196KGHOH36NDfd5P8SvUuXLuTkOB+6aZIqAm4MsYBrTBSoK2vvBx/U\nej99bFKgvLE7EV4WcI0x0UmB0sbuRHhZwDXGRCcbUjDGmAixgGuMMRESgwHXUuwECSXFTsvWF6Eh\npNgJJS1LKGWhodP3OE9tc1/uWBYNza674Hn0pSHT4IT6eUucG/U4O/cS7+a0pdgJ/c2p6co8R5lr\nYKzUmmInWtgVbj1ouZeR6jy4vCdjQ7pZPtT0PaGkk5mhv3Fc9+Myu0EnPjjty+MyO+S0Q6GkwQl1\nEovTc/+ejHVcrwkSY1e4FnCNMdHJB5Q0difCywKuMSY6xeAYrgVcY0x0soBrjDERYgHXGGMiKMYC\nrq0WZoyJTmFeLUxERovI5yKSLyKP1PB6ooisDLy+WUS6Bfa3E5H/FpHTIrIg6D3fF5Edgff8l9SR\nscCucENUXl7O/v378fmc37NpTEFBAS6XXd+ExAf8KzxViYgbWAj8GCgAtopIjqruqlJsEvCNqqaK\nyARgLnAz/nslZgDfCWxV/QGYDGzGny9tNLDmnP2wiQ/VnWviw4YNG/jXv/6FiHD55Zfj8/m4sv/3\nOXPytOO6Jd7teKKEK96NL4RJFaGUd8W78JU7/4XhjnfhdVh+9KwBrJ25xXHdofSlIT+TUM5NqOVb\nXNSKbVvzuPvuu8nLyyMtLY327duzdu3ayjJr165l2rRpeL1e7r77bh55pPoFWGlpKbfffjsff/wx\n7dq1Y+XKlXTr1s1xfxtJ/SY+dEpXpjic+PBY7RMfRGQwkKmqowLPHwVQ1aeqlFkXKLMxkAL9ENCh\nIlW6iNwJpKvq1MDzS4H/VtVvBZ5PBIaq6j3n6odd4Tpw6tQpiouLSUxMJC4ujpSUFDweD399K5sB\nAwZEpA/79+8nISGBSy65JCLtbdmy5fyP7f+FVrywsBBVJSUl5fzaC1G9ji1Ehw4doqysjMsuu4yV\nK1cyZsyY88rau2TJEtq0aUN+fj4rVqzg4YcfjvmsvUAoY7jtRaTqB7tYtdqsl07AgSrPC4CBQXVU\nllFVj4icANoBRedos1Ognqp1np2qowoLuHUoLy9nx44dJCUl4XK5UFWKi4v5+uuvSUhIoLS04deP\nU1WOHTvGZZddFpH2fD4fqhqRtiraKysri1h7IsKxY8do0aJFg7fldrs5dOgQbdq04a233qKgoIBv\nvvmGNm3aVJapmrUXqMzaWzXgZmdnk5mZCcC4ceOYOnUqqhrbSU5Du0uhyKb2NmE+n4+SkhJ8Ph9D\nhgxhy5YtqCoul4vNmzejqjRr1ox//OMfDd6X8vJyvF4vZWVlDd4W+K+4ysvLI3Js4P+sS0tLOXr0\naETa83q9/OMf/yApKSkiAcvj8ZCXl8crr7xCq1at6NSpEz169CAlJYW1a9fWmLV38+bN1eqoWiYu\nLo7k5GSOHj1K+/btG7z/jSa8t4UVAp2rPE8J7KupTEFgSCEZqO2HsjBQT211VmOj+DXw+Xxs3boV\nEal2ZetyuTh9+jQiQvPmzSPyn7Xi6i+SSQO9Xq+jrMXh4nK5IvolpNvtJj4+PmJX1HFxcTRv3pyZ\nM2fStm1bOnXqxKFDh8jPz2f06NER6UOTVDG118lWt61ALxHpLiIJwAQgOEdRDnBH4PE4YL3W8iWX\nqn4NnBSRQYG7E24Hal1gwwJukBMnTnDmzBlSU1Mrg5zH4/81e+rUKeLj4yMW/FSVkpISmjVrFtE/\nHSMdcCHyQTchIQGfz1d5bhuay+UiKSmJzMxM5s+fT7du3Th58iQ7duzgmWeeqTNrb9XMvh6PhxMn\nTtCuXbuI9L1Rhem2MFX1AFOBdcBnwCpV3Skis0SkImPnEqCdiOQDDwCV31yKyF5gHnCniBSISMV4\nz33Ay0A+8AW13KEANqRwlsTERJKSkmjbti2qitvtZtOmTZWBL5KBqKysDLfbHfHg5/P5In4Lk9vt\nxuPxkJCQELE2mzVrxr/+9S/cbndEfqFV/GVUUlLCrFmz+O1vf8vBgwfZvXs3Z86c4T/+4z94//33\nWbFiBa+//nq192ZkZPDqq68yePBg3njjDYYPHx7b47cQ9plmqroa/61bVffNrPK4BLjpHO/tdo79\neZx9q9g52RVukGbNmuFyuSqv8kpKSigpKaF58+YRDXxerxev1xvRoQRonGAL/oDr9Tq/LSscXC4X\niYmJlJREbkmqqsNRv/3tb1m2bBldu3alZcuWbNq0idatWzN+/PjKrL0VmXknTZrE0aNHSU1NZd68\necyZMydifW40FUkknWxNhN2HG+TMmTNs2bIFr9db+WVOpP+kB/99l/Hx8REPfl6vF1UlLi7yf/yU\nlpZG/BdMRbuN8Vl7PJ7KX6oPPfQQ+/fv59JLLyUxMfGs+3SbqPrdh9s+XclweB/uK7YAeZPUvHlz\nfvCDH1Q+j/lbb0yjqvj52rp1K0eOHOHiiy+2n7cKtnhN7BORiI+ZGgNEbFJLk6GEbWpvtLCAa4yJ\nTgpEdli/wVnANcZEJxtSMMaYCInBgGu3hRkTJdauXcsVV1xBampqjbd9lZaWcvPNN5OamsrAgQPZ\nu3dv5DsZSTF4W5gFXGOiQMWKYWvWrOG5557jt7/9LV26dKkWeCtWDLvvvvsoLCykb9++jBgxgn37\n9jVizxuY1+HWRFjANSYKVKwY1rVrV+6//35+9atfMXnyZLKysti1y79GdnZ2NnfccQdXXnklu3bt\nIi4ujp/85CdMnz69kXvfQMK7lkJUsIBrTBSoWA2sIvACzJs3j4KCAh544IFqZYYNG8ZFF11EcnIy\nhw8fZtWqVWetsRsTYnBIwb40MyaKFBYW0qlTJ5YtW8Z1111Hhw4deOGFF+jSpctZq5v5fD5efPFF\n2rZty/jx42nevDnf+973zlqHocmKwdvC7ArXmChQdTWwI0eOkJyczLe+9S3+/Oc/06dPH4YOHcqR\nI0dIT09nzpw5eDweCgsLad68OaWlpTz55JN07NiRHTt20KdPH1avXl1Hi01EGJNIRgO7wjUmCvTv\n3589e/YgIuzbt4+ioiK6detGcnIyqampZGdnM2zYME6ePElWVlblgvjXXnst77//PqtWrWL8+PH8\n4he/YNeuXYwZM6bp38UQg7eFWcA1JgrExcWxYMECpk2bxpdffklaWhput5sDBw7Qs2dPkpOTKS0t\nZdu2bQB88sknvPbaa/zlL3/hzJkz/O1vfyMvL4+1a9fy8ccf880337Bt2zb69evXyEdWDxVfmsUQ\nWy3MmCjz7LPPMnPmTFq2bEnnzp25/PLLyc3NZciQIbhcLt5++23KysqIj4/H4/GgqnTo0AHwD0dc\ndNFFvPDCC/zhD384K1VPhNVvtbCEdOUSh18GHmgaq4XZGK4xUWbatGlcfPHFLFq0iDZt2rBhwwYG\nDx5Mv3796Nq1Kx06dGDKlCmUlZXx85//nCuuuIKf/vSn9OzZk7S0NObNm8fTTz/N8ePH+frrrxv7\ncOonxsZwLeAaE2UqhhceeeQRcnNzueqqqzh58iTPP/88LVq0oLy8vDIFzw033MDx48fJzs4mPj6e\nI0eOcMMNN1BSUkLHjh0pLKw1p2F0i8HbwizgGhOFxowZw549e8jOziYvL4/c3FxuvfVWJk6cyNGj\nR2nVqhUAo0aN4pZbbuHw4cNs2bKF++67j8OHD1NSUhLRdEUNouK2sDDNNBOR0SLyuYjki8gjNbye\nKCIrA69vFpFuVV57NLD/cxEZVWX/XhHZISL/KyJ1j3+oaiibMaYRvPvuu9qrVy/t3LmzduzYUVVV\nZ8yYodnZ2aqqunPnTr344ou1c+fO+r3vfU/XrVunl19+uR48eLAxux1qfKm24fq+0kqdbZBXa13g\nxp/ksQeQAPwDSAsqcx/wYuDxBGBl4HFaoHwi0D1Qjzvw2l6gvdNjsitcY5qAMWPGsHv3bv72t79V\nfkE2a9YsMjL8CWfT0tJYsmQJvXv3Zvv27ZUz0S699NLKOl588UX69u1L37596d69O8OGDWuUY3HM\nh38Bcidb3QYA+ar6paqWASuAsUFlxgKvBh6/AYwIpD8fC6xQ1VJV/Qp/ht4B53NIdluYMU3ExIkT\nyc3NpaioiJSUFH73u99RXu4fwLz33nsZM2YMq1evJjU1laSkJF555ZVq77/33nu59957KS8vZ/jw\n4ZVThqNa+GaadQIOVHleAAw8VxlV9YjICaBdYP+moPdW5LFX4D0RUeCPqrq4tk5YwDWmicjKyqr1\ndRFh4cKFddYzbdo0hg8fzvXXXx+urjUc5zeitg8aQ11cV/ALk6tVtVBEOgLvi8j/qerfzlXYAq4x\nF5ClS5eyb98+FixY0NhdCbcirf0+3EKgc5XnKYF9NZUpEJE4IBk4Wtt7VbXi38Mi8hb+oYZzBlwb\nwzXmAvHxxx/z+9//ntdeey3iKeGjwFagl4h0F5EE/F+K5QSVyQHuCDweB6xX/zdjOcCEwF0M3YFe\nwBYRaSEirQBEpAUwEvi0tk7YFa4xF4gFCxZw7Nixyi/L0tPTefnllxu5V5ERGJOdCqzDf8fCn1R1\np4jMwn+HQw6wBFguIvnAMfxBmUC5VcAu/NMspqiqV0QuBt4KpLWPA15X1bW19cOm9hpjGkr9pvZK\nP4UNDksnNYmpvXaFa4yJUhVTzWKHBVxjTJSKvfUZLeAaY6KUXeEaY0yEWMA1xpgIUZzO220qLOAa\nY6KUjeEaY0yE2JCCMcZEiF3hGmNMhNgVrjHGRIhd4RpjTIRUrEAeOyzgGmOilA0pGGNMBNmQgjHG\nRIBd4RpjTIRYwDXGmAixuxSMMSZC7C4FY4yJEBtSMMaYCIm9IYULLnWnMaapqLjCdbLVTURGi8jn\nIpIvIo/U8HqiiKwMvL5ZRLpVee3RwP7PRWSU0zqDWcA1xkSpiitcJ1vtRMQNLASuAdKAiSKSFlRs\nEvCNqqYCzwFzA+9Nw5/BtzcwGlgkIm6HdVZjAdcYE6UqvjRzstVpAJCvql+qahmwAhgbVGYs8Grg\n8RvACPHnQB8LrFDVUlX9CsgP1OekzmpCHcOtV9pjY4xx7ut1kNneYeFmIpJX5fliVV1c5Xkn4ECV\n5wXAwKA6KsuoqkdETgDtAvs3Bb23U+BxXXVWY1+aGWOikqqObuw+hJsNKRhjLgSFQOcqz1MC+2os\nIyJxQDJwtJb3OqmzGgu4xpgLwVagl4h0F5EE/F+C5QSVyQHuCDweB6xXVQ3snxC4i6E70AvY4rDO\namxIwRgT8wJjslOBdYAb+JOq7hSRWUCequYAS4DlIpIPHMMfQAmUWwXswn9LxBRV9QLUVGdt/RB/\nADfGGNPQbEjBGGMixAKuMcZEiAVcY4yJEAu4xhgTIRZwjTEmQizgGmNMhFjANcaYCPn/1cwF1S6g\nCOoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3a40c0ca90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solver.pressure('pressure')\n",
    "pp.plot_grid(solver.grid(), 'pressure', color_map = [0, 0.04], if_plot=False)\n",
    "plt.title('pressure')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pressure increases in the whole domain until the shutin of the well at $t=0.05$ s, and after this the pressure eavens out. Because the fractures are more permeable than the matrix, the fluid will tend to flow out along the fractures instead of in the matrix. \n",
    "\n",
    "The time evolution of the pressure filed looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"fig/slightly_compressible.gif\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<img src=\"fig/slightly_compressible.gif\">')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
